Όλα τα απαιτούμενα δεδομένα και αρχεία θα τα βρείτε εδώ
Εννοείται ότι πριν ξεκινήσουμε να κάνουμε το οτιδήποτε, έχουμε δημιουργήσει ένα νέο project στο R-Studio, το όνομα του οποίου - για λόγους τους οποίους θα καταλάβετε αργότερα - ΔΕΝ πρέπει να περιέχει:
1. ελληνικούς χαρακτήρες
2. κενά (αντικαταστήστε τα κενά με _)
Το ίδιο ισχύει ΚΑΙ για το file path του εν λόγω φακέλου (π.χ., όχι ‘E:/Εργαστήριο Οικολογίας/SDM’, αλλά ’E:/Ergastirio_Oikologias/SDM).
Ένα άλλο κρίσιμο σημείο είναι το εξής: θα χρειαστεί να εγκαταστήσετε τα πακέτα που πρόκειται να χρησιμοποιήσετε σε αυτό το tutorial1.
Τα κλιματικά δεδομένα μας μπορούμε να τα κατεβάσουμε από το WorldClim, ενώ τα υψομετρικά δεδομένα από εδώ και τα δεδομένα σχετικά με την ξηρασία και την εξατμισιοδιαπνοή από εδώ. - Εάν φυσικά δουλεύουμε για μια συγκεκριμένη εργασία (πτυχιακή, μεταπτυχιακή, κτλ) και όχι στα πλαίσια ενός μαθήματος επιλογής. Προς το παρόν, θα χρησιμοποιήσουμε κάποιες εντολές για να κατεβάσουμε τα δεδομένα που θέλουμε μέσω της R.
Όπως μπορείτε να διαπιστώσετε, υπάρχουν πάρα πολλά κλιματικά μοντέλα (επί παραδείγματι, CCSM4, HadGEM-2, κτλ) και τέσσερα (4) διαφορετικά κλιματικά σενάρια που έχουν να κάνουν με τη συγκέντρωση των αερίων του θερμοκηπίου στην ατμόσφαιρα εντός του 21ου αιώνα { ενδεικτική βιβλιογραφία ή άλλη ενδεικτική βιβλιογραφία ή ακόμα μια ενδεικτική βιβλιογραφία ή τελευταία ενδεικτική βιβλιογραφία}2.
Συνεπώς, δεν μπορούμε να χρησιμοποιήσουμε μόνο ένα κλιματικό μοντέλο, πόσω μάλλον ένα κλιματικό σενάριο. Ευτυχώς όμως, δεν χρειάζεται και ούτε είναι επιστημονικά σωστό, να χρησιμοποιήσουμε όλα τα κλιματικά μοντέλα. Αναλόγως της περιοχής στην οποία εργαζόμαστε, κάποια κλιματικά μοντέλα αναπαριστούν καλύτερα τις περιβαλλοντικές συνθήκες σε σχέση με κάποια άλλα (McSweeney et al., 2015). Καθώς θα διαπιστώσετε (όταν με το καλό διαβάσετε την προαναφερθείσα εργασία), ακόμα και έτσι, είναι πολλά τα κλιματικά μοντέλα για την περιοχή μας. Όμως, μόνο 3 από αυτά έχουν δεδομένα και για τα 4 κλιματικά σενάρια (BCC, HadGEM, CCSM4). Εμείς χάριν συντομίας, θα εργαστούμε όμως μόνο με ένα (1) κλιματικό μοντέλο και ενα (1) κλιματικό σενάριο.
Ας φορτώσουμε τις απαραίτητες βιβλιοθήκες3.
library(easypackages)
libraries("sp", "rgdal", "raster", "dismo", "rgbif", "rgeos", "scales", "scrubr",
"mapr", "tidyverse", "rasterVis")Τώρα, ας δημιουργήσουμε τέσσερα αντικείμενα τα οποία θα περιέχουν τα αρχεία για τις παρούσες και τις μελλοντικές κλιματικές συνθήκες, τα υψομετρικά δεδομένα και την χώρα που μας ενδιαφέρει.
Πρώτα θα κατεβάσουμε τις γεωγραφικές μεταβλητές για την χώρα που μας ενδιαφέρει (οι ΗΠΑ στην προκειμένη περίπτωση). Θα πρέπει να εισάγουμε (στην εντολή country στον κώδικα που παρατίθεται πιο κάτω) τον 3-ψήφιο ISO κωδικό για τις ΗΠΑ που τον βρίσκουμε από αυτόν τον ιστότοπο ή με την εντολή ISO_countries <- getData("ISO3") %>% as.data.frame4. Εάν θέλουμε να περιορίσουμε την περιοχή μελέτης σε μια μικρότερη περιοχή, τότε θα πρέπει να χρησιμοποιήσουμε ένα άλλο λογισμικό, το Qgis, ώστε να “κόψουμε” την περιοχή που μας ενδιαφέρει. Προς το παρόν, δεν θα χρειαστεί να το κάνουμε αυτό, αλλά θα περιοριστούμε να κατεβάσουμε τα πλήρη δεδομένα για τις ΗΠΑ:
USA <- getData("GADM", country = "USA", level = 1)Το αρχείο που φτιάξαμε περιέχει κάποια (αρκετά) δεδομένα. Ας τα δούμε και ας αναπαραστήσουμε γραφικά κάποια από αυτά:
datatable(USA@data)Hawaii <- subset(USA, NAME_1 == "Hawaii")
SW_US <- subset(USA, NAME_1 == "California" | NAME_1 == "Nevada" | NAME_1 ==
"Utah" | NAME_1 == "Arizona")
plot(USA)plot(Hawaii)plot(SW_US)Κάντε το ίδιο για όποια πολιτεία της Αμερικής θέλετε. Στη συνέχεια, φτιάξτε ένα αντικείμενο το οποίο να περιέχει την Montana, το Idaho και το Wyoming και τέλος ένα αντικείμενο το οποίο να μην περιέχει την Αλάσκα και την Hawaii και ονομάστε το
USA_ter.
Μετά θα κατεβάσουμε τις παρούσες και τις μελλοντικές κλιματικές συνθήκες5, όπως αυτές δίνονται από μια παγκόσμια διαδικτυακή βάση δεδομένων, το WorldClim. Οι κλιματικές αυτές μεταβλητές είναι 19 και είναι οι ακόλουθες:
Κλιματικές μεταβλητές
| Var. | Description |
|---|---|
| BIO1 | Annual Mean Temperature |
| BIO2 | Mean Diurnal Range (Mean of monthly (max temp – min temp)) |
| BIO3 | Isothermality (BIO2/BIO7) (* 100) |
| BIO4 | Temperature Seasonality (standard deviation *100) |
| BIO5 | Max Temperature of Warmest Month |
| BIO6 | Min Temperature of Coldest Month |
| BIO7 | Temperature Annual Range (BIO5-BIO6) |
| BIO8 | Mean Temperature of Wettest Quarter |
| BIO9 | Mean Temperature of Driest Quarter |
| BIO10 | Mean Temperature of Warmest Quarter |
| BIO11 | Mean Temperature of Coldest Quarter |
| BIO12 | Annual Precipitation |
| BIO13 | Precipitation of Wettest Month |
| BIO14 | Precipitation of Driest Month |
| BIO15 | Precipitation Seasonality (Coefficient of Variation) |
| BIO16 | Precipitation of Wettest Quarter |
| BIO17 | Precipitation of Driest Quarter |
| BIO18 | Precipitation of Warmest Quarter |
| BIO19 | Precipitation of Coldest Quarter |
Όσον αφορά τις μελλοντικές κλιματικές συνθήκες, πρέπει να προσδιορίσουμε ποιό κλιματικό μοντέλο μας ενδιαφέρει (υπάρχουν πολλές επιλογές όπως προαναφέραμε, αλλά μας ενδιαφέρουν συγκεκριμένα κλιματικά μοντέλα), ποιό κλιματικό σενάριο (υπάρχουν 4 διαφορετικά: 26, 45, 60, 85), καθώς και ποιά χρονολογία μας ενδιαφέρει (50 ή 70 - αναφέρεται στα έτη 2050 και 2070, αντίστοιχα). Προς το παρόν εμείς θα δουλέψουμε με το μοντέλο CCMS4 (αντιστοιχεί στο μοντέλο CC), το κλιματικό σενάριο 85 και το έτος 2070. Τις μελλοντικές κλιματικές συνθήκες θα τις χρησιμοποιήσουμε προκειμένου να δείτε πώς θα προβλέψουμε την κατανομή των ειδών στον χώρο σε άλλο χρόνο:
current_climate <- getData("worldclim", var = "bio", res = 2.5)
current_climate## class : RasterStack
## dimensions : 3600, 8640, 31104000, 19 (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667 (x, y)
## extent : -180, 180, -60, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## names : bio1, bio2, bio3, bio4, bio5, bio6, bio7, bio8, bio9, bio10, bio11, bio12, bio13, bio14, bio15, ...
## min values : -278, 9, 8, 64, -86, -559, 53, -278, -501, -127, -506, 0, 0, 0, 0, ...
## max values : 319, 213, 96, 22704, 489, 258, 725, 376, 365, 382, 289, 10577, 2437, 697, 265, ...
future_climate <- getData("CMIP5", var = "bio", res = 2.5, rcp = 85, model = "CC",
year = 70)
future_climate## class : RasterStack
## dimensions : 3600, 8640, 31104000, 19 (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667 (x, y)
## extent : -180, 180, -60, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs
## names : cc85bi701, cc85bi702, cc85bi703, cc85bi704, cc85bi705, cc85bi706, cc85bi707, cc85bi708, cc85bi709, cc85bi7010, cc85bi7011, cc85bi7012, cc85bi7013, cc85bi7014, cc85bi7015, ...
## min values : -217, 9, 7, 54, -58, -491, 53, -284, -436, -94, -443, 0, 0, 0, 0, ...
## max values : 351, 219, 96, 21959, 536, 281, 707, 426, 411, 430, 325, 13161, 3017, 682, 223, ...
Ας αλλάξουμε τα ονόματα των μεταβλητών μας ως ακολούθως:
names(current_climate) <- c("bio_1", "bio_2", "bio_3", "bio_4", "bio_5", "bio_6",
"bio_7", "bio_8", "bio_9", "bio_10", "bio_11", "bio_12", "bio_13", "bio_14",
"bio_15", "bio_16", "bio_17", "bio_18", "bio_19")
names(future_climate) <- c("bio_1", "bio_2", "bio_3", "bio_4", "bio_5", "bio_6",
"bio_7", "bio_8", "bio_9", "bio_10", "bio_11", "bio_12", "bio_13", "bio_14",
"bio_15", "bio_16", "bio_17", "bio_18", "bio_19")Καθώς τα κλιματικά μας δεδομένα όπως διαπιστώσατε έχουν εξαιρετικά μεγάλο εύρος τιμών, θα χρειαστεί να τα μετατρέψουμε σε oC όσον αφορά την θερμοκρασία (bio_1). Εδώ μπορείτε να διαπιστώσετε το γιατί6. Αυτό γίνεται με την ακόλουθη εντολή:
gain(current_climate$bio_1) = 0.1
gain(future_climate$bio_1) = 0.1Στη συνέχεια, θα κατεβάσουμε τα υψομετρικά δεδομένα για την χώρα που μας ενδιαφέρει:
altitude <- getData("alt", country = "USA", mask = TRUE)
altitude## $`C:/R_SDM/Command_Lines/SDM/hsdm-master/SDM/Correct files for SDM LAB EMB/USA1_msk_alt.grd`
## class : RasterLayer
## dimensions : 3012, 6948, 20927376 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -124.8, -66.9, 24.4, 49.5 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84
## data source : C:\R_SDM\Command_Lines\SDM\hsdm-master\SDM\Correct files for SDM LAB EMB\USA1_msk_alt.grd
## names : USA1_msk_alt
## values : -171, 4275 (min, max)
##
##
## $`C:/R_SDM/Command_Lines/SDM/hsdm-master/SDM/Correct files for SDM LAB EMB/USA2_msk_alt.grd`
## class : RasterLayer
## dimensions : 2448, 5916, 14482368 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -179.2, -129.9, 51.1, 71.5 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84
## data source : C:\R_SDM\Command_Lines\SDM\hsdm-master\SDM\Correct files for SDM LAB EMB\USA2_msk_alt.grd
## names : USA2_msk_alt
## values : -22, 6098 (min, max)
##
##
## $`C:/R_SDM/Command_Lines/SDM/hsdm-master/SDM/Correct files for SDM LAB EMB/USA3_msk_alt.grd`
## class : RasterLayer
## dimensions : 228, 900, 205200 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : 172.4, 179.9, 51.2, 53.1 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84
## data source : C:\R_SDM\Command_Lines\SDM\hsdm-master\SDM\Correct files for SDM LAB EMB\USA3_msk_alt.grd
## names : USA3_msk_alt
## values : -1, 1099 (min, max)
##
##
## $`C:/R_SDM/Command_Lines/SDM/hsdm-master/SDM/Correct files for SDM LAB EMB/USA4_msk_alt.grd`
## class : RasterLayer
## dimensions : 420, 672, 282240 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -160.3, -154.7, 18.8, 22.3 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84
## data source : C:\R_SDM\Command_Lines\SDM\hsdm-master\SDM\Correct files for SDM LAB EMB\USA4_msk_alt.grd
## names : USA4_msk_alt
## values : -3, 4120 (min, max)
## We are only interested in the first of the four elements of the
## altitudinal data
altitude <- altitude[[1]]Στη συνέχεια, θα ‘κόψουμε’ τις περιβαλλοντικές μας μεταβλητές στο σχήμα των ΗΠΑ και ακολούθως θα δημιουργήσουμε ένα νέο αρχείο που θα περιέχει τις μεταβλητές αυτές και την περιοχή μελέτης μας:
cropped_current <- crop(current_climate, USA_ter, snap = "in") %>% mask(., USA_ter)
cropped_future <- crop(future_climate, USA_ter, snap = "in") %>% mask(., USA_ter)Τώρα, θα ‘κόψουμε’ τις κλιματικές μεταβλητές για τις ΗΠΑ σε επιμέρους αρχεία:
SW_crop_current <- crop(cropped_current, SW_US, snap = "in") %>% mask(., SW_US)
SW_crop_future <- crop(cropped_future, SW_US, snap = "in") %>% mask(., SW_US)Κάντε το ίδιο για την California και το Oregon.
Μετά, μπορούμε να αναπαραστήσουμε γραφικά τα αποτελέσματα:
gplot(cropped_current[[1]]) + geom_raster(aes(fill = value)) + scale_fill_gradientn(colours = c("brown",
"red", "yellow", "darkgreen", "green")) + coord_equal() + xlab("Longitude") +
ylab("Latitude") + ggtitle("Temperature")Κάντε το ίδιο για το αρχείο που φτιάξατε προηγουμένως.
Στη συνέχεια, θα υπολογίσουμε κάποια βασικά στατιστικά στοιχεία για κάθε αρχείο το οποίο δημιουργήσαμε:
cellStats(cropped_current, "quantile")## bio_1 bio_2 bio_3 bio_4 bio_5 bio_6 bio_7 bio_8 bio_9 bio_10 bio_11
## 0% -7.1 59 22 1794 60 -245 120 -135 -157 5 -157
## 25% 6.6 123 31 7431 278 -130 350 105 -44 184 -56
## 50% 9.8 134 37 8385 305 -85 382 175 36 213 -13
## 75% 14.6 155 41 9626 329 -28 419 211 161 248 46
## 100% 25.2 213 65 13060 456 178 509 332 329 360 212
## bio_12 bio_13 bio_14 bio_15 bio_16 bio_17 bio_18 bio_19
## 0% 46 7 0 5 19 0 1 11
## 25% 392 66 11 20 170 39 121 51
## 50% 692 96 21 36 262 77 212 112
## 75% 1067 121 53 56 331 187 296 233
## 100% 3345 556 154 104 1512 480 660 1455
cellStats(cropped_current, "mean")## bio_1 bio_2 bio_3 bio_4 bio_5 bio_6
## 10.59007 138.48284 36.36840 8517.96950 301.83563 -79.58974
## bio_7 bio_8 bio_9 bio_10 bio_11 bio_12
## 381.42537 153.05698 55.71655 213.60417 -5.79786 747.70502
## bio_13 bio_14 bio_15 bio_16 bio_17 bio_18
## 98.39876 32.30871 38.61828 267.33779 113.08449 209.89919
## bio_19
## 161.71277
cellStats(cropped_current, "sd")## bio_1 bio_2 bio_3 bio_4 bio_5 bio_6
## 5.268230 21.448783 6.533989 1729.940968 38.814827 69.275880
## bio_7 bio_8 bio_9 bio_10 bio_11 bio_12
## 54.956739 83.319657 113.217157 44.299732 67.841672 415.811072
## bio_13 bio_14 bio_15 bio_16 bio_17 bio_18
## 51.187839 26.819703 20.246943 144.910739 88.540570 112.976747
## bio_19
## 151.904632
Κάντε το ίδιο για το SW_crop και για το αρχείο που φτιάξατε για την California και το Oregon.
Τέλος, μπορούμε να δημιουργήσουμε υποσύνολα των περιοχών που μας ενδιαφέρουν, ώστε να ανταποκρίνονται σε συγκεκριμένες κλιματικές (ή άλλες συνθήκες):
hot_US <- cropped_current[[1]] < 25 & cropped_current[[1]] > 14.6
plot(hot_US)Δοκιμάστε το και για άλλες περιοχές και μεταβλητές.
Τώρα, θα πρέπει να ελέγξουμε εάν τα υψομετρικά δεδομένα έχουν την ίδια ανάλυση (resolution), μέγεθος (nrow, ncol) και έκταση (bbox), με τα περιβαλλοντικά μας δεδομένα:
altitude## class : RasterLayer
## dimensions : 3012, 6948, 20927376 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -124.8, -66.9, 24.4, 49.5 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84
## data source : C:\R_SDM\Command_Lines\SDM\hsdm-master\SDM\Correct files for SDM LAB EMB\USA1_msk_alt.grd
## names : USA1_msk_alt
## values : -171, 4275 (min, max)
cropped_current## class : RasterBrick
## dimensions : 596, 1387, 826652, 19 (nrow, ncol, ncell, nlayers)
## resolution : 0.04166667, 0.04166667 (x, y)
## extent : -124.75, -66.95833, 24.54167, 49.375 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : bio_1, bio_2, bio_3, bio_4, bio_5, bio_6, bio_7, bio_8, bio_9, bio_10, bio_11, bio_12, bio_13, bio_14, bio_15, ...
## min values : -7.1, 59.0, 22.0, 1794.0, 60.0, -245.0, 120.0, -135.0, -157.0, 5.0, -157.0, 46.0, 7.0, 0.0, 5.0, ...
## max values : 25.2, 213.0, 65.0, 13060.0, 456.0, 178.0, 509.0, 332.0, 329.0, 360.0, 212.0, 3345.0, 556.0, 154.0, 104.0, ...
Στην περίπτωση κατά την οποία τα δύο σετ δεδομένων έχουν διαφορετική ανάλυση, έκταση και μέγεθος, θα χρειαστεί να αλλάξουμε την ανάλυση των υψομετρικών δεδομένων7.
agg_altitude <- aggregate(altitude,
5, ## Why choose 5 and not any other number?
fun = mean) Μετά, θα αλλάξουμε το μέγεθος των pixel των υψομετρικών δεδομένων για να συμφωνεί με αυτό των περιβαλλοντικών δεδομένων:
elev_US_resample <- resample(agg_altitude, cropped_current, method = "bilinear")Τέλος, θα τα ενώσουμε σε ένα αρχείο, τόσο για το παρόν, όσο και για το μέλλον:
var_US <- readAll(stack(cropped_current, elev_US_resample))
var_US_70 <- readAll(stack(cropped_future, elev_US_resample))
elev_SW_US <- crop(elev_US_resample, SW_US, snap = "in") %>% mask(., SW_US)
var_SW_US <- readAll(stack(SW_crop_current, elev_SW_US))
var_SW_70 <- readAll(stack(SW_crop_future, elev_SW_US))Καλό είναι να αποθηκεύσετε τα αρχεία αυτά υπό την μορφή Rds και να τα αφαιρέσετε από την μνήμη της R8. Αυτό το κάνουμε, προκειμένου να μπορούμε εύκολα να το φορτώνουμε στην R και να μην υπάρχουν τα δεδομένα αυτά συνεχώς στην μνήμη της session που έχουμε ανοικτή, διότι αυξάνουν αρκετά τον χρόνο απόκρισης της (καθυστερεί πολύ η R). Μπορείτε ανά πάσα στιγμή να τα φορτώσετε στην R, με την ακόλουθη εντολή:
## Creating a folder to save our data
dir.create("./Data/Variables", recursive = TRUE, showWarnings = FALSE)
## Save our data in rds format
saveRDS(var_US, "./Data/Variables/Environmental variables and Elevation (30s) USA.rds")
saveRDS(var_US_70, "./Data/Variables/Environmental variables and Elevation (30s) USA 2070.rds")
## Load our data from rds format
var_US <- readRDS("./Data/Variables/Environmental variables and Elevation (30s) USA.rds")
var_US_70 <- readRDS("./Data/Variables/Environmental variables and Elevation (30s) USA 2070.rds")Ας αναπαραστήσουμε γραφικά το αποτέλεσμα, το οποίο θα φανεί στην ακόλουθη εικόνα:
plot(var_SW_US[[c(1, 12, 20)]], col = rainbow(100, start = 0, end = 0.8))Τρέξτε την εντολή
cellStatsγια τα 3 αρχεία τα οποία έχετε φτιάξει (var_US, var_SW_US και var_CO_US). Βλέπετε διαφορές;
Ας αναπαραστήσουμε γραφικά - αλλά με τέτοιον τρόπο ώστε να μας δίνει κάποιες στατιστικές πληροφορίες - τις περιβαλλοντικές μας μεταβλητές:
histogram(var_SW_US[[20]])levelplot(var_SW_US[[20]])Κάντε το ίδιο και για άλλες μεταβλητές
Ας δούμε τα δεδομένα μας τρισδιάστατα:
plot3D(var_SW_US[[20]])Κάντε το ίδιο και για άλλες μεταβλητές
Σε αυτήν την ενότητα θα διαπιστώσουμε εάν οι περιβαλλοντικές μας μεταβλητές έχουν υψηλή συσχέτιση μεταξύ τους στην περιοχή μελέτης.
library(usdm)Στο σημείο αυτό θα μετατρέψουμε τα raster αρχεία μας σε data-frames, προκειμένου να μπορέσουμε να διαπιστώσουμε ποιές μεταβλητές εμφανίζουν υψηλή συσχέτιση:
var_US_df <- na.omit(as.data.frame(var_US))
var_SW_US_df <- na.omit(as.data.frame(var_SW_US))Γιατί το κάνουμε αυτό μόνο για τις παρούσες κλιματικές συνθήκες;
Όπως θα γνωρίζετε, ΔΕΝ μπορούμε να χρησιμοποιήσουμε σε μια ανάλυση δύο μεταβλητές, οι οποίες εμφανίζουν συντελεστή συσχέτισης \(> 0.7\)9, καθότι στην προκειμένη περίπτωση οι δύο αυτές μεταβλητές είναι συγγραμικές και ως εκ τούτου εάν τις συμπεριλάβουμε και τις δύο στην ανάλυση μας, τότε δεν θα γνωρίζουμε σε ποιά εκ των δύο θα οφείλεται το αποτέλεσμα της ανάλυσης μας.
Προκειμένου να ξεπεράσουμε αυτόν τον σκόπελο, θα τρέξουμε την ακόλουθη εντολή10:
vifcor(var_SW_US_df, th = 0.7)## 13 variables from the 20 input variables have collinearity problem:
##
## bio_16 bio_6 bio_19 bio_12 bio_10 bio_17 bio_11 bio_7 bio_3 bio_5 bio_1 bio_14 bio_15
##
## After excluding the collinear variables, the linear correlation coefficients ranges between:
## min correlation ( bio_18 ~ bio_13 ): 0.02776225
## max correlation ( bio_13 ~ bio_4 ): -0.6901746
##
## ---------- VIFs of the remained variables --------
## Variables VIF
## 1 bio_2 1.281558
## 2 bio_4 2.403638
## 3 bio_8 2.160976
## 4 bio_9 1.807053
## 5 bio_13 2.232166
## 6 bio_18 2.482672
## 7 USA1_msk_alt 4.043632
Τώρα δοκιμάστε το και εσείς, αλλάζοντας το όριο σε 0.8 και 0.9. Τι παρατηρείτε; Ποιές μεταβλητές θα κρατούσατε στην ανάλυση;
Πλεόν είμαστε σε θέση να επιλέξουμε τις μη συγγραμικές μεταβλητές, τις οποίες θα χρησιμοποιήσουμε στην ανάλυση μας.
## If these are our uncorrelated variables
var_SW_US_sub <- stack(subset(var_SW_US, c(2, 3, 8, 13:15, 20)))
var_SW_70_sub <- stack(subset(var_SW_70, c(2, 3, 8, 13:15, 20)))
## Let's scale our variables (mean = 0, sd = 1)
var_SW_US_sub <- stack(scale(var_SW_US_sub))
var_SW_70_sub <- stack(scale(var_SW_70_sub))Κάντε το ίδιο για το σύνολο των ΗΠΑ.
Στο σημείο αυτό, είμαστε πλέον στην ευχάριστη θέση να μπορούμε να χρησιμοποιήσουμε τα δεδομένα θέσης για τα είδη τα οποία μας ενδιαφέρουν. Υπάρχουν δύο περιπτώσεις, όσον αφορά την απόκτηση τέτοιων δεδομένων:
- να τα έχουμε αποκτήσει πρωτογενώς (δηλαδή να έχουμε εργαστεί στο πεδίο και να έχουμε εντοπίσεις τις θέσεις εμφάνισης του εκάστοτε είδους και των υποπληθυσμών του)
- να τα έχουμε αποκτήσει δευτερογενώς (δηλαδή είτε να έχουμε χρησιμοποιήσει τις διαδικτυακές βάσεις ψηφιοποιημένων - και ενίοτε γεωαναφερμένων - δεδομένων GBIF και BIEN, είτε να έχουμε ψηφιοποιήσει χάρτες εμφάνισης θέσης11)
Το ιδανικό βεβαίως είναι να έχουμε αποκτήσει πρωτογενώς τα στοιχεία αυτά, καθότι στην περίπτωση αυτή, μπορούμε να είμαστε σίγουροι για τις γεωγραφικές συντεταγμένες των υποπληθυσμών του προς μελέτη είδους. Στην πλειονότητα των περιπτώσεων όμως, δεν έχουμε πρόσβαση σε τέτοια δεδομένα, οπότε αναγκαστικά καταφεύγουμε στην δεύτερη λύση και δη στις διαδικτυακές βάσεις δεδομένων. Ευτυχώς για εμάς, υπάρχουν κάποια πακέτα12 στην R, όπως το rgbif και το rbien, τα οποία αυτοματοποιούν την διαδικασία μεταφόρτωσης δεδομένων από τις διαδικτυακές αυτές βάσεις13.
Εμείς θα εργαστούμε με το πακέτο rgbif προς το παρόν.
library(easypackages)
libraries("viridis", "ggmap", "biogeo", "spThin", "mapr")Σήμερα θα ασχοληθούμε με ένα είδος το οποίο είναι ενδημικό της δυτικής ακτής των ΗΠΑ, το Yucca brevifolia.
Το είδος Yucca brevifolia στην πολιτεία της Nevada
Σύμφωνα με το Υπουργείο Αγροτικής Ανάπτυξης των ΗΠΑ, το είδος αυτό είναι ενδημικό των πολιτειών της California, της Nevada, της Utah και της Arizona. Ανήκει στην οικογένεια Agavaceae και το κοινό του όνομα είναι Joshua tree, εξαιτίας της μορφής του: οι Μορμόνοι άποικοι οι οποίοι αντίκρυσαν για πρώτη φορά αυτό το δένδρο στην έρημο Mojave τον 18ο αιώνα θέωρησαν ότι τους θύμιζε τον ικέτη Ιησού που σηκώνει τα χέρια του προς τον ουρανό καθώς προσεύχεται.
Ένα άλλο κοινό όνομα του είδους αυτού είναι το izote de desierto (Ισπανικά: το μαχαίρι της ερήμου). Η Yucca brevifolia εντοπίζεται κυρίως στην έρημο Mojave σε υψόμετρο μεταξύ 400 έως και 1800 m, είναι ένα ταχυαυξές είδος (7.6 cm/έτος για τα πρώτα 10 χρόνια) το οποίο δεν σχηματίζει ετήσιους δακτύλιους. Το ρίζικο του σύστημα φθάνει σε βάθος 11 m. Τα ψηλότερα άτομα έχουν 15 m ύψος, ενώ η περίοδος ανθοφορίας του είναι μεταξύ Φεβρουαρίου και Απριλίου: όπως τα περισσότερα είδη τα οποία απαντώνται σε ερήμους, χρειάζονται τη σωστή ποσότητα βροχής το κατάλληλο χρονικό διάστημα, προκειμένου να ανθίσουν. Έχουν δε, βρεθεί άτομα ηλικίας 1000 ετών. Σύμφωνα με τους Shafer et al. (2001),η κλιματική αλλαγή αναμένεται να έχει αρνητική επίδραση στο είδος αυτό.
Πρώτα χρειάζεται να βεβαιωθούμε ότι το προς μελέτη είδος, υπάρχει όντως εντός της διαδικτυακής βάσης την οποία επιθυμούμε να χρησιμοποιήσουμε, προκειμένου να εξάγουμε δεδομένα. Η εντολή name_suggest() ερευνά όλα τα είδη (taxa καλύτερα) τα οποία μοιάζουν με το είδος το οποίο ψάχνουμε και κρατά μόνο το όνομα των taxa το οποίο ταιριάζει ακριβώς με αυτό το οποίο μας ενδιαφέρει14.
spp_int <- name_suggest(q = "Yucca brevifolia", rank = "species", limit = 10000)
(spp_int <- spp_int[grepl("^Yucca brevifolia$", spp_int$canonicalName), ])| key | canonicalName | rank |
|---|---|---|
| 7520449 | Yucca brevifolia | SPECIES |
| 7636266 | Yucca brevifolia | SPECIES |
| 2775592 | Yucca brevifolia | SPECIES |
| 7390928 | Yucca brevifolia | SPECIES |
Στο βήμα αυτό αναγνωρίσαμε τον μοναδικό αριθμό αναφοράς (key) του προς μελέτη είδους και μόνο χρησιμοποιώντας τον αριθμό αυτό μπορούμε να είμαστε βέβαιοι ότι τα δεδομένα θέσης τα οποία θα λάβουμε από την διαδικτυακή αυτή βάση δεδομένων θα αναφέρονται στο είδος το οποίο μας ενδιαφέρει.
sp_occ <- occ_search(taxonKey = spp_int$key, country = "US", fields = c("name",
"key", "country", "decimalLatitude", "decimalLongitude"), hasCoordinate = T,
limit = 1000, return = "data")
sp_occ## $`7520449`
## [1] "no data found, try a different search"
##
## $`7636266`
## [1] "no data found, try a different search"
##
## $`2775592`
## # A tibble: 1,000 x 5
## name key decimalLongitude decimalLatitude country
## <chr> <int> <dbl> <dbl> <chr>
## 1 Yucca brevifolia 1806362560 -116 33.9 United St~
## 2 Yucca brevifolia 1807289573 -116 34.1 United St~
## 3 Yucca brevifolia 1807298279 -114 36.6 United St~
## 4 Yucca brevifolia 1807303045 -116 34.1 United St~
## 5 Yucca brevifolia 1806341722 -115 36.2 United St~
## 6 Yucca brevifolia 1805394434 -116 34.1 United St~
## 7 Yucca brevifolia 1831095189 -116 34.0 United St~
## 8 Yucca brevifolia 1805399185 -116 34.0 United St~
## 9 Yucca brevifolia 1805417949 -114 36.7 United St~
## 10 Yucca brevifolia 1831185811 -114 36.5 United St~
## # ... with 990 more rows
##
## $`7390928`
## [1] "no data found, try a different search"
##
## attr(,"args")
## attr(,"args")$hasCoordinate
## [1] TRUE
##
## attr(,"args")$limit
## [1] 1000
##
## attr(,"args")$offset
## [1] 0
##
## attr(,"args")$taxonKey
## [1] 7520449 7636266 2775592 7390928
##
## attr(,"args")$country
## [1] "US"
##
## attr(,"args")$fields
## [1] "name" "key" "country"
## [4] "decimalLatitude" "decimalLongitude"
Προκειμένου να μην προκύψουν προβλήματα με το μονοπάτι (path) του φακέλου και των δεδομένων μας, καλό είναι να αφαιρέσουμε τα κενά εντός του ονόματος του είδους μας και να διαλέξουμε το ID15 του taxon μας με τις περισσότερες εγγραφές:
sp_occ <- sp_occ[["2775592"]]
sp_occ$name <- sub(" ", ".", sp_occ$name)
(spp_to_model <- unique(sp_occ$name))## [1] "Yucca.brevifolia"
Τώρα έχει έρθει η στιγμή να δούμε πόσα γεωαναφερμένα δεδομένα θέσης έχουμε στην διάθεση μας:
sort(table(sp_occ$name), decreasing = T)## Yucca.brevifolia
## 1000
Ας οπτικοποιήσουμε τα αποτελέσματα της ανάλυσης μας έως αυτό το σημείο:
## Setting the extent (data taken from our RasterStack)
myloc <- c(-125, -67, 24.6, 49.3)
mymap <- get_map(location = myloc,
source = 'stamen',
maptype = 'terrain',
crop = F)
## Vizualise our map
ggmap(mymap) +
geom_point(aes(x = decimalLongitude, y = decimalLatitude),
data = sp_occ,
alpha = .5,
color = "darkred",
size = 3) +
xlab('Longitude') +
ylab('Latitude') Θα πρέπει να τονίσουμε το εξής: Εφ’όσον δεν έχουμε συλλέξει πρωτογενώς τα δεδομένα μας, θα χρειαστεί να ελέγξουμε εάν οι θέσεις εμφάνισης που μεταφορτώσαμε από τις διαδικτυακές βάσεις δεδομένων, συμφωνούν με τα γνωστά όρια εξάπλωσης του είδους που μας ενδιαφέρει. Εάν όντως υπάρχουν κάποια σημεία τα οποία βρίσκονται εκτός της περιοχής εξάπλωσης, θα χρειαστεί να τα αφαιρέσουμε από την ανάλυση. Θα χρησιμοποιήσουμε το πακέτο scrubr για να κάνουμε αυτή την ανάλυση16.
## Let's create a new object, containing only the correct coords
sp_occ_ck <- dframe(sp_occ) %>% coord_impossible() %>% coord_incomplete() %>%
coord_unlikely()
## Is there any difference between the two sets of coords?
NROW(sp_occ)## [1] 1000
NROW(sp_occ_ck)## [1] 1000
Ας φτιάξουμε έναν φάκελο ώστε να αποθηκεύσουμε τα δεδομένα μας και ας τα σώσουμε, προκειμένου να μπορέσουμε να τα φορτώσουμε κάποια άλλη στιγμή:
dir.create("./Species", recursive = TRUE, showWarnings = FALSE)
dir.create("./Species/Yucca brevifolia", recursive = TRUE, showWarnings = FALSE)
saveRDS(sp_occ_ck, "./Species/Yucca brevifolia.rds")Πλέον μπορούμε να συνδέσουμε τις γεωαναφερμένες θέσεις εμφάνισης του είδους που μας ενδιαφέρει με τα raster αρχεία τα οποία είχαμε δημιουργήσει προηγουμένως. Στο σημείο αυτό, μπορούμε να εξάγουμε τα περιβαλλοντικά δεδομένα για κάθε μια εκ των θέσεων εμφάνισης του είδους που μας ενδιαφέρει με μια μόνο εντολή, αναδεικνύοντας κατ’αυτόν τον τρόπο το πλεονέκτημα του να έχουμε ενώσει προηγουμένως τα raster αρχεία μας σε ένα RasterStack.
coordinates(sp_occ_ck) <- c("longitude", "latitude")
pts.clim <- raster::extract(var_SW_US, sp_occ_ck, method = "bilinear")
sp_occ_clim <- data.frame(cbind(coordinates(sp_occ_ck), pts.clim, sp_occ_ck@data))Ακολούθως, θέλουμε να αφαιρέσουμε τις θέσεις εμφάνισης οι οποίες πέφτουν εκτός των περιβαλλοντικών μας μεταβλητών (μπορεί να πέφτουν στην θάλασσα).
sp_occ_clim <- sp_occ_clim %>% na.omit()
coordinates(sp_occ_clim) <- c("longitude", "latitude")Και τώρα, ήρθε η στιγμή να οπτικοποιήσουμε τα αποτελέσματα μας και δη, την παρουσία του είδους μας στην περιοχή μελέτης:
map.ext <- extent(var_SW_US)
plot(var_SW_US[[20]], col = terrain.colors(100), ext = map.ext)
plot(sp_occ_clim, pch = 16, cex = 0.5, add = T)
plot(SW_US, add = T)Στο σημείο αυτό, θα ήθελα να σας τονίσω ότι ορισμένες φορές, πέρα από τους ελέγχους που διενεργήσαμε έως τώρα, ίσως χρειαστεί να προχωρήσουμε σε μια αραίωση (thinning) των θέσεων εμφάνισης. Αυτό μπορεί να είναι απαραίτητο όταν έχουμε πάρα πολλά σημεία τα οποία είναι αρκετά κοντά το ένα στο άλλο17. Σε κάθε περίπτωση, το επιστημονικά ορθό είναι να υπάρχει μια (1) θέση εμφάνισης ανά km2. Ένα πακέτο που μπορεί να μας βοηθήσει είναι το spThin18.
## Min. distance: 1000 m Number of reps: 1000
spoccs <- correct_occurrence_points_sp %>% as.data.frame() %>% rename(x = coords.x1,
y = coords.x2, Species = name) %>% addmainfields(species = "Species")
thinned <- thin(loc.data = spoccs, lat.col = "y", long.col = "x", spec.col = "Species",
thin.par = 1, reps = 1000, verbose = T, out.dir = "./Thinned", locs.thinned.list.return = TRUE,
write.files = TRUE, out.base = "Yucca brevifolia thinned")McSweeney, C.F., Jones, R.G., Lee, R.W., & Rowell, D.P. (2015) Selecting CMIP5 GCMs for downscaling over multiple regions. Clim. Dyn., 44, 3237–3260.
Shafer, S.L., Bartlein, P.J., & Thompson, R.S. (2001) Potential changes in the distributions of western north america tree and shrub taxa under future climate scenarios. Ecosystems, 4, 200–215.
Υπάρχουν δύο τρόποι για να το κάνετε αυτό γρήγορα:
1. library(easypackages) packages('package_name_1', 'package_name_2')
2. install.packages(c('package_name_1', 'package_name_2'), dependecies = TRUE)↩
Καλό είναι να τις διαβάσετε↩
Εάν δεν τις έχετε εγκαταστήσει, κάντε το τώρα: install.packages(insert_package_name_here) ή library(easypackages), packages('raster', 'sp', 'rgdal') - εάν θέλετε να εγκαταστήσετε τα τρια αυτά πακέτα↩
Ή με αυτήν την εντολή εάν θέλουμε μόνο τις ΗΠΑ: ISO_countries <- getData("ISO3") %>% as.data.frame %>% filter(NAME=="USA")↩
Επειδή θα πάρει αρκετό χρόνο με τη σύνδεση στο εργαστήριο, να τις έχετε κατεβάσει. Διαφορετικά, ζητήστε μου να σας τις δώσω εγώ↩
Γιατί;↩
Γιατί όμως να μην αλλάξουμε την ανάλυση των περιβαλλοντικών δεδομένων;↩
Πώς τα αφαιρούμε από την μνήμη της R;↩
Αυτό είναι το λεγόμενο hard boundary - άλλοι ερευνητές εφαρμόζουν το όριο \(> 0.8\) ή ακόμα και \(> 0.9\), αλλά είναι επιστημονικά ασφαλέστερο να εφαρμόζει κανείς το αυστηρότερο όριο↩
Εδώ εφαρμόζουμε το αυστηρότερο όριο - μπορείτε να το αλλάξετε, εάν επιθυμείτε↩
Ένα τυπικό παράδειγμα είναι ο Άτλαντας της Χλωρίδας του Αιγαίου ή ο Άτλαντας της Χλωρίδας της Ευρώπης, όμως στην περίπτωση αυτή υπεισέρχονται άλλα προβλήματα τα οποία σχετίζονται με την ανάλυση των εν λόγω χαρτών (π.χ., ο Άτλαντας της Χλωρίδας της Ευρώπης είναι σε ανάλυση \(50 \times 50 km^2\))↩
Δεν είναι τα μοναδικά όμως: υπάρχουν πολλά άλλα, όπως το dismo και το spocc↩
Κάθε ένα από αυτά τα πακέτα, έχει τουλάχιστον μια σύντομη περιγραφή (αγγλιστί: vignette), την οποία είναι ΑΠΑΡΑΙΤΗΤΟ να διαβάσετε - ειδικά αυτές του rgbif↩
Με την εντολή ?rgbif::name_suggest() μπορείτε να δείτε τα arguments που δέχεται η εντολή αυτή και τι μπορείτε να αλλάξετε - σας συνιστώ εντόνως να το κάνετε, αλλά και για κάθε εντολή που χρησιμοποιείτε↩
Εάν υπάρχουν πολλά βέβαια↩
Δείτε το manual του πακέτου scrubr ή τρέξτε τις εντολές μια μια για να δείτε τι κάνουν στον κώδικα που έχετε στα χέρια σας↩
Σε απόσταση \(<1000m\) όταν εργαζόμαστε με χάρτες υποβάθρου κλίμακας 1 x 1 km2 (ή ανάλυσης 30 arc-seconds)↩